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ABSTRACT 

This paper describes the development of an aeroelastic code (TURBO-AE) based on 
an Euler / Navier-Stokes unsteady aerodynamic analysis. A brief review of the 
relevant research in the area of propulsion aeroelasticity is presented. The paper 
briefly describes the original Euler / Navier-Stokes code (TURBO) and then details 
the development of the aeroelastic extensions. The aeroelastic formulation is 
described. The modeling of the dynamics of the blade using a modal approach is 
detailed, along with the grid deformation approach used to model the elastic 
deformation of the blade. The work-per-cycle approach used to evaluate aeroelastic 
stability is described. Representative results used to verify the code are presented. 
The paper concludes with an evaluation of the development thus far, and some 
plans for further development and validation of the TURBO-AE code. 


INTRODUCTION 

NASA’s Advanced Subsonic Technology (AST) program seeks to develop 
technologies to increase the fuel efficiency of commercial aircraft engines, improve 
the safety of engine operation, reduce the emissions, and reduce engine noise. With 
the development of new designs of ducted fans, compressors, and turbines to achieve 
these goals, a basic aeromechanical requirement is that there should be no flutter or 
high resonant blade stresses in the operating regime. In order to verify the 
aeroelastic soundness of the design, an accurate prediction of the unsteady 
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aerodynamics and structural dynamics of the propulsion component is required. The 
complex geometry, the possibility of shock waves and flow separation makes the 
modeling of the unsteady aerodynamics a difficult task. The advanced blade 
geometry, new blade materials and new blade attachment concepts, make the 
modeling of the structural dynamics a difficult problem. 

Computational aeroelastic modeling of fans, compressors, and turbines, requires 
several simplifying assumptions. Flutter calculations are typically carried out 
assuming that the blade row is isolated. This simplifies the structural dynamics 
formulation and the unsteady aerodynamic calculations considerably. 

For an isolated blade row flutter calculation, the modeling of the unsteady 
aerodynamics is the biggest challenge. Many simplifying assumptions are made in 
the modeling of the unsteady aerodynamics. In the past, a panel method based on 
linear compressible small-disturbance potential theory has been used to model the 
unsteady aerodynamics and aeroelasticity of propfans in subsonic flow [Williams, 
1990, and Kaza et al., 1989]. The major limitation of this analysis is the neglect of 
transonic and viscous flow effects in the model. Although this analysis requires a 
fairly small computational effort, the inherent limitation in the model precludes its 
use in a majority of practical applications. 

More recently, a full potential unsteady aerodynamic analysis has been used with a 
modal structural dynamics approach to model the aeroelastic behavior of fan blades 
[Bakhle and Reddy, 1993 and Bakhle et al., 1993]. In this aeroelastic analysis, the 
unsteady aerodynamics model is in the time domain. The time domain flutter 
calculations use the method of simultaneous integration of structural dynamics and 
aerodynamics equations in time. Alternatively, for the frequency domain flutter 
calculations, a Fourier analysis is required to transform the time domain unsteady 
aerodynamic force coefficients to the frequency domain. The eigenvalue approach is 
then used. Although the full potential aerodynamic formulation is able to model 
transonic effects to a certain extent (limited to weak shocks), the vortical effects are 
still neglected. Thus, for example, the blade tip vortex, or a leading-edge vortex is 
not modeled. Significant computational effort is required for such flutter 
calculations based on full potential aerodynamics. An aeroelastic analysis based on 
the Euler equations [Srivastava and Reddy, 1995] can model vortical flows, but 
viscous effects are not taken into account. Recently, other researchers [Hall, 1992, 
He and Denton, 1993, Gerolymos and Vallet, 1994, Peitsch et al., 1994, and 
Carstens, 1994] have also developed inviscid and viscous unsteady aerodynamic 
analyses for vibrating blades. 
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For aeroelastic problems in which viscous effects play an important role (such as 
flutter with flow separation, or stall flutter, and flutter in the presence of shock and 
boundary-layer interaction), a more advanced aeroelastic computational capability 
is required. This paper describes the development of such an aeroelastic code 
(TURBO-AE). This aeroelastic code is based on an unsteady aerodynamic Euler / 
Navier-Stokes code (TURBO), developed separately, which is described in the 
following section. In this paper, the grid deformation method is described, along 
with the interpolation from the finite-element structural mesh to the CFD mesh, 
and the calculation of the aeroelastic forces and work. 


DESCRIPTION OF EULER / NAVIER-STOKES CODE- TURBO 

This section describes i>erv briefly the TURBO code. Additional details regarding the 
code are available elsewhere [Janus, 1989 and Chen, 1991]. The TURBO code 
provides all the unsteady aerodynamics to the TURBO-AE code. 

The TURBO code was originally developed [Janus, 1989] as an in viscid flow solver 
for modeling the flow through multistage turbomachinery. It has the capability to 
handle multiple blade rows with even or uneven blade count, stationary or rotating 
blade rows and blade rows at an angle of attack. Multiple blade passages are 
included in the calculation, when required. Additional developments were made 
[Chen, 1991] to incorporate viscous terms into the model. The code can now be 
applied to model realistic turbomachinery configurations with flow phenomena such 
as shocks, vortices, separated flow, secondary flows, and shock and boundary layer 
interactions. 

The code is based on a finite volume scheme. Flux vector splitting is used to 
evaluate the flux Jacobians on the left hand side of the governing equations [Janus, 
1989] and Roe’s flux difference splitting is used to form a higher-order TVD (Total 
Variation Diminsihing) scheme to evaluate the fluxes on the right hand side. 
Newton sub-iterations are used at each time step to maintain higher accuracy. A 
Baldwin-Lomax algebraic turbulence model is used in the code. 


DEVELOPMENT OF AEROELASTIC CODE - TURBO-AE 

The TURBO-AE code assumes a normal mode representation of the structural 
dynamics of the blade. Thus, the dynamic characteristics of each blade are assumed 
to be represented in terms of in-vacuum modes, with the associated natural 
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frequency and generalized mass for each mode. Typically, a finite-element analysis 
code such as NASTRAN is used to calculate the modal data mentioned. No 
restriction is placed on the analysis, except that it should provide grid point 
coordinates, modal deflections, natural frequencies, and generalized mass for each 
mode of interest. 

At the current stage of development, the aeroelastic code (TURBO-AE) models only 
in-phase blade vibrations. That is, the interblade phase angle is restricted to be 
zero. This is a limitation of the aeroelastic code and not of the original unsteady 
aerodynamics code (TURBO). Thus, in TURBO-AE, only one blade and one blade 
passage are modeled. 

A work-per-cycle approach is used to determine aeroelastic (flutter) stability. Using 
this approach, the motion of the blade is prescribed to be a harmonic vibration in a 
specified in-vacuum normal mode with a specified frequency. The vibration 
frequency is typically the natural frequency for the mode of interest, but some other 
frequency can also be used. The aerodynamic forces acting on the vibrating blade 
and the work done by these forces on the vibrating blade during a cycle of vibration 
are calculated. If work is being done on the blade by the aerodynamic forces, the 
blade is dynamically unstable, since it will result in extraction of energy from the 
flow, leading to an increase in amplitude of oscillation of the blade. Note that 
coupled mode flutter cannot be modeled with this approach. 

In the following two sub-sections, the grid deformation scheme and the work-per- 
cycle and generalized force calculations are detailed. 


Grid Deformation 

The blade of a rotor undergoes a rigid-body rotation about the axis of the rotor and a 
simultaneous elastic deformation. This is modeled through grid point motion. This 
requires the calculation of a new grid at each time step. The grid point motion due 
to blade rotation consists simply of a rotation of each grid point about the axis of the 
rotor. There is no relative motion between grid points. The grid point motion due to 
the elastic deformation of the blade is superposed on this rotation. The elastic 
deformation of the blade does cause relative motion between grid points and 
therefore, the term grid deformation is used in this context. 

The grid deformation approach [Huff, 1989] used in the TURBO-AE code is 
described here. In this approach, the grid points on the surface of the blade move by 
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an amount equal to the deformation of the blade. The grid points on the remaining 
boundaries of the computational domain do not move. The motion of the grid points 
in the interior of the computational domain is determined as a product of three 
coefficients (one for each computational coordinate direction) and the motion of 
points on blade surfaces. The coefficient in each computational direction is linearly 
dependent on the distances of the interior grid point from the boundaries and from 
the blade surfaces along that computational direction. Thus, points close to the 
moving blade move almost as much as points on the blade, and points close to the 
stationary boundaries do not move much at all. The coefficients are not presented 
here, but can be found in [Huff, 1989], where they are referred to as weighting 
functions. 

The grid deformation is calculated as described here. First, new locations of all grid 
points on the blade surface are calculated. For a harmonic vibration of the blade in a 
selected in-vacuum normal mode, the displacement of any point on the blade (due to 
elastic deformation) X(x,y,z,t) can be written in terms^of the generalized coordinate 
q(t) and the modal deflection or mode shape function d(x,y^z,t) as: 


X(x,y,z,t) = q(t) h(x,yjz,t) (1) 

Note that x, y, and z represent the coordinates in a rotating coordinate system. 
Further note that the modal deflection or mode shape function 8 has been 
interpolated from the finite-element grid to the CFD grid. This interpolation, which 
is described in the following paragraph, is performed only once for each mode of 
interest. 

The interpolation of the modal deflections from the finite element grid onto the CFD 
grid is accomplished as follows. The interpolation is done with the blade undeflected 
in a reference position. At each CFD grid point on the blade surface, the distance to 
the nearest three finite-element grid points is calculated. Then, the modal 
deflections at these three nearest neighbors are used in a bi-linear interpolation 
scheme to calculate the interpolated valu^of the modal deflection at that CFD grid 
point. The interpolated modal deflections 8 are stored and used at each time step to 
calculate the motion of each grid point on the blade surface. 

Next, the new locations of the interior points are calculated. This requires the 
calculation of coefficients based on the reference location of the interior grid point. 
As previously mentioned, the distance of each interior grid point, from the 
computational boundaries and from the blade surfaces, along the corresponding 
computational direction, is calculated. For example, for an interior grid point (i,j, k), 
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the distance to the nearest i = constant (inlet/exit) boundary is calculated along the 
{-direction coordinate line. Then, the distance from the interior point to the relevant 
blade surface ({-index corresponding to leading/trailing edge) is calculated along the 
{-direction coordinate line. These two calculations are repeated for the other two 
computational directions. Then, these distances are used to calculate the coefficients 
in the three computational coordinate directions. The motion of the interior grid 
points is simply a linear scaling of the motion of the blade surface. 

Special attention is required for grid points located between the blade tip and the 
duct or casing. Typically, the grid points in this region are very closely spaced in the 
radial direction. If it is assumed that the grid points on the casing do not move when 
the blade deforms, the computational cells can become excessively skewed due to 
the motion of the blade tip and the interior grid points. To avoid this, the grid points 
on the casing are allowed to move (along the casing surface) in such a way as to 
reduce the skewing of the interior computational cells. The location of the grid 
points in the interior region is calculated by linear interpolation between the blade 
tip and the casing. This treatment improves the aspect ratio of the computational 
cells considerably. 


Work and Force Calculation 

To determine aeroelastic stability using the work-per-cycle approach, the blade 
motion is specified to be harmonic: 


q(t) =qosin(<£t) 

where qo is the amplitude of motion and co is the vibration frequency. 
The work-per-cycle done on the blade is calculated as: 


W 


= <j> j p dA • (dx/dt)dt 


or, 


M 


W = <b |p<£A*5go® cos(ast) dt 

s 


( 2 ) 


( 3 ) 


( 4 ) 
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In the above, p = p(x,y,z,t ) is the unsteady pressure on the blade surface due to 
blade vibration, A is the blade surface area vector pointing into the blade surface, 
j is the integral over the blade surface, j> is the integral over one cycle of blade 
S 

vibration. For viscous calculations, the viscous stresses must be included in the 
calculation along with the pressure. This is done in the TURBO-AE code, but the 
expressions for the work including inviscid and viscous contributions are not 
presented here. 

The work-per-cycle is an indicator of aeroelastic stability. The blade is dynamically 
unstable if the work done qr the blade during a cycle of blade vibration is positive. 

Finally, the generalized force on the blade A m>n is defined as: 

A m ,n — j Pn dA • 5/jj (5) 

S 

where p n is the unsteady pressure due to blade vibration in the vfi 1 mode and b m is 
the 171 th modal deflection. The viscous terms are included in the generalized force 
calculation, along with the pressure term. However, the viscous contributions are 
not presented here for simplicity. The generalized force can be used to determine 
the flutter stability of the blade. For an analysis with just one mode, the blade will 
flutter if the generalized force for that mode leads the motion. This provides a 
simple check of the work-per-cycle calculation. 


SAMPLE RESULTS 

In this section, some sample results are presented. The configuration selected is 
derived from the Energy Efficient Engine (E-cubed) fan rotor. The E-cubed program 
was established by GE Aircraft Engines under NASA sponsorship in the 1980’s to 
demonstrate component technologies necessary to achieve higher efficiencies and 
reduce environmental effects in future subsonic turbofan engines. A summary of the 
program has been presented recently [Smith, 1993]. Details regarding design and 
performance tests have also been presented [Sullivan and Hager, 1983]. 

The results presented here are meant to demonstrate the state of development of 
the code. The fan rotor has 32 blades with a tip diameter of 210.8 cm (83 inches). 
The inlet flow (axial) Mach number used in this calculation is 0.5. The CFD grid is a 
simple grid with 15 points on the blade surface in both the chordwise and spanwise 
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directions. It is assumed that the tip gap is zero. The finite-element structural 
dynamics data is for a grid with 224 points. The results presented are for an in viscid 
run of the code. The code has been run in the viscous mode for other configurations. 

To begin, a steady solution is obtained for this configuration. The aeroelastic 
calculations are then performed. Figure 1 shows the finite-element and CFD grids 
on the blade surface. It is to be noted that the blade planforms from the two grids 
are slightly different. In the regions where the two planforms do not match, 
extrapolation (rather than interpolation) is used to transfer the modal deflections 
from the finite-element grid to the CFD grid, and hence there is potential for errors 
at these locations 

Figure 2 shows the original and interpolated modal deflections corresponding to the 
first mode. The quantity plotted is the magnitude of the deflection at each grid 
point. Note that the original data on the finite-element grid is provided on a mean 
surface, whereas the interpolated data is for the CFD grid on the two distinct blade 
surfaces. Overall, the interpolated deflections match the original deflections. 
However, some differences are noted in the blade tip region. Referring to Figure 1, it 
can be noted that this is one of the areas in which the blade planforms from the 
CFD grid and the finite-element grid do not match. Figure 3 shows the original and 
interpolated modal deflections for the second mode. In this case, no si gnifican t 
disagreement can be seen between the original and interpolated data. 

Aeroelastic calculations have been performed for the first two modes separately. 
These two modes have natural frequencies of about 72 Hz and 164 Hz respectively. 
Aeroelastic calculations were performed for these modes at their natural vibration 
frequencies, using 100 computational time steps per cycle of blade vibration. Note 
that, the time step used in the two calculations is not the same. In Figure 4, the 
instantaneous work is plotted against the time for blade vibration in the first mode. 
The amplitude of blade vibration used in this calculation results in a maximum 
deflection of the blade of about 0.9% of the blade tip diameter. The code was run for 
eight cycles of blade vibration. The variation of instantaneous work with time shows 
that the flow has become periodic in time. The variation is seen to be nearly 
sinusoidal, indicating the absence of significant non-linear effects (as evidenced by 
the absence of higher harmonic content). This may be attributed to the subsonic 
flowfield. Also, the results indicate that the selected amplitude of blade vibration 
does not result in any non-linear effects. 

In Figure 5, the variation of the cumulative work, from the beginning of the current 
vibration cycle, is plotted. This quantity starts at zero at the beginning of each 
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vibration cycle of the blade. At the end of the cycle, the cumulative work done is the 
work-per-cycle. This is represented as a symbol at the end of each cycle of vibration. 
After eight cycles, it can be seen that the work-per-cycle remains negative. Thus, 
work is being done tx the blade, and this shows that the blade is stable under the 
conditions of the analysis. 

Figure 6 shows the work-per-cycle after each cycle of vibration. This information is 
also seen in Figure 5, although on a slightly different scale. It is presented to show 
the convergence of the flow to periodicity in time. For the configuration analyzed, it 
can be seen that the work-per-cycle does not vary much after the fourth cycle of 
blade vibration. 

Figure 7 shows the variation of the generalized force on the blade (An). The motion 
of the blade, given by Equation (2), is also shown. It may be seen that the 
generalized force An lags the motion. For an analysis with just one mode, this 
indicates that the blade vibrations are stable, since a force lagging the motion will 
not result in work being done qr the blade. This is a verification of the conclusion 
reached from the work-per-cycle calculation. 

Finally, results are presented for the blade vibrating in the second mode. 
Calculations are performed for an amplitude of blade vibration that results in a 
mavirmim deflection of the blade of about 0.2% blade tip diameter. Figure 8 shows 
the instantaneous work plotted against the time. Although the instantaneous work 
becomes periodic in time after about four or five cycles, the variation is seen to be 
significantly non-linear (as evidenced by the higher harmonic content). This is in 
contrast to the nearly linear variation obtained for the first mode (Figure 4) for a 
comparatively larger maximum blade deflection. The cause of this non-linear 
behavior is unknown and is being investigated. The response is expected to be linear 
for a sufficiently small amplitude. This is seen in Figure 9, where the amplitude of 
vibration is reduced in half (0.1% blade tip diameter). It should be noted that the 
work-per-cycle approach is valid for either amplitude, since the approach only 
requires a periodic response, not necessarily a linear one. 

Figure 10 shows the work-per-cycle for blade vibration in the second mode for the 
smaller amplitude. It can be seen that the work-per-cycle converges to a small 
negative value. Note that the work-per-cycle is not normalized by the amplitude of 
vibration. Hence, the smaller values work-per-cycle are partly a result of the 
smaller amplitude. The calculation indicates that the blade vibrations in the second 
mode are stable. It should be noted that the work-per-cycle continues to change 
slightly, even after eight cycles of blade vibration. 
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CONCLUDING REMARKS 


An aeroelastic analysis code named TURBO-AE has been developed. The starting 
point for the development was an Euler / Navier-Stokes unsteady aerodynamic code 
named TURBO. Routines have been developed to interpolate the structural 
deflections from the finite-element grid to the CFD grid. Grid deformation routines 
have been developed to calculate a new grid for the deformed blade at each time 
step. Routines have been developed for the calculation of work and generalized 
forces. These routines have been verified by running the code for a realistic 
configuration. 

Results have been presented to show the working of the code for two modes of blade 
vibration. Dynamically, the blade is seen to be stable in both the first and second 
modes. The results for the first mode show that the force acting on the blade is 
linear at the selected amplitude of vibration. The results for the second mode show 
that a non-linear response is obtained for a fairly small amplitude of vibration. 
However, a reduction of amplitude results in a linear response, as expected. 

The calculations for the second mode illustrate an advantage of the work-per-cycle 
approach, namely, it remains valid even if the force is non-linear. However, the 
work-per-cycle approach suffers from an assumption that the aerodynamics and the 
structural dynamics can be decoupled. 

A major limitation of the code is that it is currently restricted to the analysis of in- 
phase blade motions. In a propulsion component, fan, compressor, or turbine, it is 
necessary to consider the various interblade phase angles at which flutter can occur. 
Hence, in the future, the code will be extended to allow the analysis of arbitrary 
interblade phase angles. This can be accomplished either by including multiple 
blade passages in the calculations, or by using a single blade passage with time (or 
phase) shifted boundary conditions. Also, it is necessary that the TURBO-AE code 
be exercised to evaluate its ability to analyze and predict flutter for conditions in 
which viscous effects are significant. This is also planned for the future. 
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Figure 1 : Typical finite-element and CFD grids, and superimposed planforms. 
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Figure 2: Original and interpolated modal deflections for first mode; interpolated 

data is for two blade surfaces. 
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Figure 3: Original and interpolated modal deflections for second mode; 
interpolated data is for two blade surfaces. 
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Figure 4: Instantaneous work done on blade vibrating in first mode. 
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Figure 5: Cumulative work (from beginning of each cycle) done on blade vibrating in 
first mode. Symbols indicate the value at the end of each cycle. 
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Figure 6: Work-per-cycle for blade vibrating in first mode. 
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Figure 7: Generalized force acting on blade vibrating in first mode. 



Figure 8: Instantaneous work done on blade vibrating in second mode with 
maximum blade deflection 0.2% blade tip diameter. 
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Figure 9: Instantaneous work done on blade vibrating in second mode with 
maximum blade deflection 0.1% blade tip diameter. 



Figure 10: Work-per-cycle for blade vibrating in second mode. 


17 



REPORT DOCUMENTATION PAGE 


Form Approved 
OMB No. 0704-0188 


Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, 
gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this 
collection of information, including suggestions for reducing this burden, to Washington Headquarters Services, Directorate for Information Operations and Reports, 1215 Jefferson 
Davis Highway, Suite 1204, Arlington, VA 22202-4302, and to the Office of Management and Budget, Paperwork Reduction Project (0704-0188), Washington, DC 20503 


1. AGENCY USE ONLY {Leave blank) 


4. TITLE AND SUBTITLE 


2. REPORT DATE 


November 1996 


3. REPORT TYPE AND DATES COVERED 

Technical Memorandum 


Development of an Aeroelastic Code Based on an Euler/Navier- Stokes 
Aerodynamic Solver 


6. AUTHOR(S) 

Milind A.Bakhle, Rakesh Srivastava, Theo G. Keith, Jr, George L. Stefko, and 
J. Mark Janus 


5. FUNDING NUMBERS 


WU-538-06-14 

NAG3-1803 


7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 

National Aeronautics and Space Administration 
Lewis Research Center 
Cleveland, Ohio 44135-3191 


8. PERFORMING ORGANIZATION 
REPORT NUMBER 


E- 10523 


9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 

National Aeronautics and Space Administration 
Washington, DC 20546-0001 


10. SPONSORING/MONITORING 
AGENCY REPORT NUMBER 


NASA TM- 107362 


11. SUPPLEMENTARY NOTES 

Milind A. Bakhle, Rakesh Srivastava, and Theo G. Keith, Jr., University of Toledo, Toledo, Ohio; George L. Stefko, 
NASA Lewis Research Center, Cleveland, Ohio; J. Mark Janus, Mississippi State University, Mississippi State, Missis- 
sippi. Responsible person, George L. Stefko, organization code 5230, (216) 433-3920. 


12a. DISTRIBUTION/AVAILABILITY STATEMENT 


12b. DISTRIBUTION CODE 


Unclassified - Unlimited 
Subject Category 39 

This publication is available from the NASA Center for AeroSpace Information, (301) 621-0390. 


13. ABSTRACT (Maximum 200 words) 

This paper describes the development of an aeroelastic code (TURBO-AE) based on an Euler/Navier- Stokes unsteady 
aerodynamic analysis. A brief review of the relevant research in the area of propulsion aeroelasticity is presented. The 
paper briefly describes the original Euler/Navier-Stokes code (TURBO) and then details the development of the 
aeroelastic extensions. The aeroelastic formulation is described. The modeling of the dynamics of the blade using a modal 
approach is detailed, along with the grid deformation approach used to model the elastic deformation of the blade. The 
work-per-cycle approach used to evaluate aeroelastic stability is described. Representative results used to verify the code 
are presented. The paper concludes with an evaluation of the development thus far, and some plans for further develop- 
ment and validation of the TURBO-AE code. 


14. SUBJECT TERMS 


Aeroelasticity; Flutter; Viscous; Inviscid; Navier-Stokes 


17. SECURITY CLASSIFICATION 
OF REPORT 

Unclassified 


18. SECURITY CLASSIFICATION 
OF THIS PAGE 

Unclassified 


19. SECURITY CLASSIFICATION 
OF ABSTRACT 

Unclassified 


15. NUMBER OF PAGES 

19 


16. PRICE CODE 

A03 


20. LIMITATION OF ABSTRACT 


NSN 7540-01-280-5500 


Standard Form 298 (Rev. 2-89) 
Prescribed by ANSI Std. Z39-18 
298-102 



























